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Abstract 


The hydrological design of the reservoirs to be built or the revision of the 
existing ones, requires the estimation of the so-called Design Flood 
Hydrograph. The simplest and most approximate way to estimate such a 


hydrograph, for a determined joint return period, is through the bivariate 
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Frequency Analysis of the maximum flow (Q) and the annual runoff 
volume (V) of the registered floods. Copula Functions (FC) are 
probabilistic models based on the dependence between Q and V, which 
easily establish their bivariate distribution, based on previously adopted 
marginal functions or distributions of any type, equal or different. The 
application of the FC in the hydrological estimates showed that a decisive 
aspect in their ideal selection is related to the dependence on the extreme 
right of the data (AF) and that which have (A,) certain FC. Therefore, in 
this study the FC are exposed: Student's t, Gumbel-Hougaard, Clayton 
Associate and Joe, which show increasing values Ay. The values of 2, are 
contrasted against the 15% obtained in 16 random real records of Q and 
V, to establish the applicability of each cited FC. In addition, the record of 
26 annual Q and V data of the inflow floods to the Adolfo López Mateos 
Dam (Humaya), in the state of Sinaloa, Mexico, is processed as a 
numerical application. Finally, the Conclusions are presented, which 
highlight the advantages of FC in the bivariate frequency analysis of 


floods. 


Keywords: Copula Functions, Kendall's tau quotient, Spearman's rho 
coefficient, dependence at the upper tail, observed dependence, joint 


return periods, secondary return period. 


Resumen 


El diseño hidrológico de los embalses por construir o la revisión de los ya 
existentes requiere la estimación del llamado hidrograma de la creciente 


de diseño. La manera más simple y aproximada para estimar tal 
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hidrograma, para un determinado periodo de retorno conjunto, es a 
través del análisis de frecuencias bivariado del gasto máximo (Q) y el 
volumen escurrido (VW) anuales de las crecientes registradas. Las 
funciones Cópula (FC) son modelos probabilísticos basados en la 
dependencia entre Q y V, que establecen fácilmente su distribución 
bivariada con base en las funciones marginales previamente adoptadas o 
distribuciones de cualquier tipo, iguales o diferentes. La aplicación de las 
FC en las estimaciones hidrológicas mostró que un aspecto decisivo en su 
selección idónea está relacionado con la dependencia en el extremo 
derecho de los datos (A£'“) y la que tienen (Ay) ciertas FC. Por lo anterior, 
en este estudio se exponen las FC: t de Student, Gumbel-Hougaard, 
Clayton Asociada y Joe, que muestran valores de 1, que van en aumento. 
Se contrastan los valores de A, contra los 15“ obtenidos en 16 registros 
reales aleatorios de Q y V, para establecer la aplicabilidad de cada FC 
citada. Además, se procesa como aplicación numérica el registro de 26 
datos de Q y V anuales de las crecientes de entrada a la presa Adolfo 
López Mateos (Humaya) del estado de Sinaloa, México. Por último, se 
exponen las conclusiones, las cuales destacan las ventajas de las FC en 


los análisis de frecuencias bivariados de crecientes. 


Palabras clave: funciones Cópula, cociente tau de Kendall, coeficiente 
rho de Spearman, dependencia en el extremo superior, dependencia 


observada, periodos de retorno conjuntos, periodo de retorno secundario. 
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Introduction 


Generalities 


The Mexican Republic is located within the area of hurricane influence 
which originate in the Caribbean Sea and the Pacific Ocean. In addition, 
it is affected by other meteorological phenomena of great spatial scope, 
such as cold fronts and extensive convective storms, some of orographic 
origin. The aforementioned atmospheric phenomena produce rains of 
great magnitude, which generate floods or maximum flows that inundate 
extensive regions of the country. Such floods can cause loss of human 
lives and tremendous economic and environmental damages. This 
scenario that repeats every year, in different areas of Mexico, emphasizes 
the importance of studying the floods, in order to understand their 
dynamics and try to quantify their magnitudes, with the purpose of 
designing building and protecting hydraulic works for protection and 
utilization (Aldama, 2000; Aldama, Ramírez, Aparicio, Mejía-Zermeño, 8 
Ortega-Gil, 2006; Gómez, Aparicio, 8, Patiño, 2010). 


The main hydraulic protection works are the dams and flood bariers, 
which are generally related to rectifications and channeling. Bridges and 
other urban drainage works are also protection structures, whose 
hydrological design is based on the Design Floods (DF), which are 
maximum river flows associated with low probabilities of being exceeded. 


The most reliable estimation of the DFs is made through the Flood 
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Frequencies Analysis (AFC, by its acronym in Spanish), which tries to 
interpret the record of maximum annual flows available, in terms of future 
events of a certain probability of occurrence (Bobée 8, Ashkar, 1991; Rao 
8: Hamed, 2000; Meylan, Favre, € Musy, 2012). 


For the results of the AFC to be approximate and relatively reliable, 
the following four requirements must be met: (1) the recording of 
maximum annual flows must be random; (2) the probability distribution 
function (PDF) that represents it must be suitable; (3) the fit of the tested 
PDFs must be carried out with efficient statistical methods and (4) the 
selection of results must be objective (Hosking €: Wallis, 1997; Stedinger, 
2017). 


On the other hand, the hydrological design of reservoirs or storage 
dams requires the estimation of the so-called Hydrograph of the Design 
Flood, since in such hydraulic works the volume of the flood is extremely 
important to establish its dangerousness and define its management and 
temporary storage without risk. Aldama (2000) indicated that the 
reservoirs are not sensitive to the time that elapses in reaching the 
maximum flow, but to its magnitude and the volume of the flood. The 
later makes it possible to significantly simplify the multivariate analysis of 
the floods, since the bivariate approach of the annual maximum flow (Q) 
and volume (V) allows the hydrograph of the design flood to be defined 
with sufficient approximation, since V has a positive correlation with the 
total duration (D) of the hydrograph. 


Bivariate flood analysis began at the end of the last century, with the 
work of Goel, Seth and Chandra (1998), Yue (1999), and Yue, Ouarda, 


Bobée, Legendre and Bruneau (1999), who processed the data sets of Q 
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and V and of V and D. At the beginning of this century, Yue and 
Rasmussen (2002) established the basic hydrological concepts of 
bivariate probabilistic analysis. The first flood bivariate analysis used the 
bivariate Normal distribution and the Logistic model that accepts marginal 
functions of extreme values such as the Gumbel and GVE distributions. 
These models have two important constraints: (1) their marginal 
distributions are equal and (2) the adjustment parameters of such 
marginals intervene in the joint modeling of the Q and V register (Favre, 
El Adlouni, Perreault, Thiémonge, 8: Bobée, 2004). 


Currently, such constraints are avoided by applying the Copula 
Functions (CF), mathematical models that accept any PDF as marginal, to 
form the bivariate or multivariate distribution and which are based on the 
correlation that exists between Q and V (Shiau, Wang, 8 Tsai, 2006; 
Zhang € Singh, 2006). The use of CFs has been consolidated, as proven 
by specialized texts (Salvadori, De Michele, Kottegoda, € Rosso, 2007; 
Zhang € Singh, 2019; Chen 8 Guo, 2019; Chowdhary € Singh, 2019). 


Ramírez-Orozco and Aldama (2000), Escalante-Sandoval and Reyes- 
Chávez (2002), as well as Volpi and Fiori (2012), and Requena, Mediero 
and Garrote (2013) indicate that the bivariate AFC leads to an infinity of 
peak flow and volume combinations for an adopted joint exceedance 
probability. This implies that for the same joint return period there are 
many floods or hydrographs that will produce different effects on the 
reservoir under design or revision; logically deciding for the one that 
generates the most severe conditions in its storage and spillway of 


excesses. 
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Objectives 


The works of Dupuis (2007), and Poulin, Huard, Favre and Pugin (2007) 
showed that a decisive aspect in the selection of the ideal CF to use is 
related to the dependence on the right extreme of the data and the one 
shown by certain CF families. Therefore, the basic objectives of this study 
were the following four: [1] expose four families of CF that have 
significant dependence on their right tail (Ay); [2] describe how the Au is 
estimated from the available Q and V record (AF); [3] to contrast the 
values of AF“ obtained in 16 real Q and V records, against those which 
define the four exposed CF and [4] to process as a numerical application 
of 26 data records of annual Q and V of floods at the entrance to the 


Adolfo López Mateos Dam (Humaya), in the state of Sinaloa, Mexico. 


Copula Functions 


Concept and definition 


As already indicated, the key advantage of the Copula Functions (CF) 
consists in allowing to express a joint distribution of correlated random 
variables, as a function of their previously adopted marginal distributions. 
So, a CF links or relates the univariate marginal distributions to form a 
multivariate distribution. Another basic advantage of CFs when forming 


multivariate distributions is the fact that they separate the effects of 
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dependency among the random variables from those of the marginal 


distributions in joint modeling. 


Due to the above, the construction of the multivariate distribution is 
reduced to the study of the relationship between the correlated variables, 
provided the univariate marginal distributions are known. The use of CF 
offers complete freedom to adopt or select the univariate marginal 
distributions which best represent the data (Shiau et a/., 2006; Kottegoda 
8: Rosso, 2008; Meylan et a/., 2012; Zhang €: Singh, 2019). 


Since in this study CF will be applied in the bivariate frequency 
analysis of annual floods, the following definition refers to two correlated 
random variables X and Y, whose joint cumulative probability distribution 
function is Fx yv(x,y) with univariate marginal distributions Fx(x) and Fx(y); 
then the CF exists, it is C[:] and it is such that: 


FxyGco y) = CIF), Fy(y)] (1) 


The above equation defines the basic concept for the development 
of CFs and is known as Sklar's Theorem exposed in 1959 (Shiau et al., 
2006; Kottegoda 8 Rosso, 2008; Meylan et a/., 2012; Zhang € Singh, 
2019). 


Tecnología y ciencias del agua, ISSN 2007-2422, 


Open Access bajo la licencia ee BY-NC-SA 4.0 14(5), 120-188. DOI: 10.24850/j-tyca-14-05-03 
(https: //creativecommons.org/licenses/by-nc-sa/4.0/) 


a 0) Check for updates 
OPEN ACCESS 


Tecnología y 


CienciaszAgua 
Copula Functions Families 


The Copula functions (CF) that have been developed have been classified 
into four classes: Archimedean, extreme value, elliptic, and 
miscellaneous. They are also classified as single-parameter or multi- 
parameter copulas, depending on the extent to which the dependency 
structure between the variables Q and V is defined (Meylan et al., 2012; 
Chowdhary 8 Singh, 2019). 


Salvadori et al. (2007) provide a comprehensive and useful summary 
of CFs that have been applied in the field of hydrology. In this regard, in 
the following relationship the families of Ali-Mikhail-Haq Copulas, called 
AMH Copula, have not been included, since they restrict the dependence 
to a Kendall tau quotient which ranges from -0.1817 to 0.3333 and the 
Farlie-Gumbel-Morgenstern family, designated FGM, which restricts it 
further, to the interval of +0.2222. 


Archimedean copulas have found broad application due to their 
simple construction, single parameter, wide range, and acceptance of 
both types of dependency (positive and negative). Designating Fx(x)=u, 
FNy)=v and 0 the parameter that measures the dependence or 
association between u and v, we have the following two families of 
Archimedean Copulas, which accept negative and positive dependence 
(Shiau et al., 2006; Genest 8 Favre, 2007; Salvadori et a/., 2007; Michiels 
8 De Schepper, 2008; Zhang € Singh, 2019; Chen € Guo, 2019; 
Chowdhary 8 Singh, 2019). 
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1, Clayton. This Copula function is called Cook-Johnson by Zhang and 
Singh (2006). Sraj, Bezak and Brilly (2015) cite it with both names. 


Its equation and variation space of O are: 
_(,/-8 -9 -1/0 
C(u,v) = (u +v*"-— 1) [—1, 00) X f0) (2) 


For the positive dependence 6 > O and for the negative -1<0<0, 
with O = O for the independence between u and v. The relationship of O 


with the Kendall tau quotient is as follows: 
E (3) 


2. Frank. Its equation and variation space of 6 are: 


e79U-1)(e70 
e79-1 


CQu,0) == In 14 | (oo, 00) 110) (4) 


For the negative dependence 0<8<1 and for the positive O > 1, 
with O = 1 for the independence between u and v. The relationship of O 


with r,, is as follows: 


tn =1+5[D1(0)—1] (5) 
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where Di(0) is the Debye function of order 1, whose expression ¡s: 


DO) =5b) ds (6) 


es-1 


The previous equation was estimated by means of numerical 
integration, ratifying its results with the values tabulated by Stegun 
(1972). 


Next, a family of Extreme Value Copulas is cited, which accepts only 


positive dependence. 
3. Gumbel-Hougaard. Its equation and variation space of 6 are: 
C(u,v) = exp [E Inu)? + In a [1, 00) (7) 


With O = 1 we have independence between u and v. The relationship 
of O to Kendall's tau ratio is as follows: 


n= (8) 


Lastly, two families of Miscellaneous class Copulas are cited. 
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4. Plackett accepts negative (0 < 1) and positive (6 > 1) dependence, 


with O = 1 for independence between u and v: 


E — . 
(unas la, MPA ARA 


a 0-0 con0<g<o (9) 


The relationship of O to Spearman's rho is as follows: 


_(0+1)  201n(0) 
n= 0) e-1? 0) 


5. Raftery accepts only positive dependence and 6 varies between O 


and 1, with O = O for independence between u and v: 


C(u,v) = min(u, v) + =- (uv) /0=9)(1 — [max(u, 9] eBAan (11) 


The relationship of O to Kendall's tau quotient and to Spearman's 
rho coefficient are as follows: 


— (12) 
_— 0(4-30) 
Pn 7 (2-09)? (13) 


Tecnología y ciencias del agua, ISSN 2007-2422, 


Open Access bajo la licencia eel BY-NC-SA 4.0 14(5), 120-188. DOI: 10.24850/j-tyca-14-05-03 
(https://creativecommons.org/licenses/by-nc-sa/4.0/) 


a 0) Check for updates 
OPEN ACCESS 


Tecnología y 


CienciaszAgua 


The families of Elliptic Copulas are distributions and their marginals 
that accept the property of being elliptic, such as the Normal, Student's t 
and Cauchy distributions. A random vector X e R* is said to have an elliptic 
distribution if it admits the following stochastic representation: X =u+ 
RAU; in which, ye R*, R is an independent random variable of U, which is 
a uniform random vector on the unit sphere in R“ and A is a fixed dxd 
matrix, such that A-A' is nonsingular (Salvadori et a/., 2007; Chen 8 Guo, 
2019). 


When the marginal distributions are of different types and possibly 
not elliptic, the elliptic Copulas are called Meta-elliptic Copulas, whose 
main families are Gaussian and Student's t (Chowdhary 8 Singh, 2019). 
The Student's t-Copula will be described in detail later. 


Families of associated Copulas 


The four existing classes of Copula families and in the case of the bivariate 
ones, have three associated families, defined by the following 
transformations (Salvadori 8. De Michele, 2004; Michiels 8 De Schepper, 
2008; Chowdhary 8 Singh, 2019): 


C'(uv) =u-C(u,1—v) (14) 


C"(uv)=w-C(1-—u,v) (15) 
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The first two transformations change the nature of the dependency 
from negative to positive and vice versa. The third transformation is called 
the Copula of Survival (Survival Copula), because it is related to the joint 
distribution of the exceedance probabilities, of the pair of random 
variables u and v whose Copula is C. Furthermore, all three families of 
associated Copulas have the same value of Kendall's tau ratio; the first 


two with their sign changed. 


Dupuis (2007) highlighted the importance of the Associated 
Clayton Copula, as it has a significant dependence on its right tail defined 
by the following Equation (17) and whose formula, according to equations 
(2) and (16), is Equation (18): 


Ay =2 9 (17) 


Cu,v) =u+v-1+[U1-0 + (1-w) 9-1] (18) 


Numerical integration 


To quantify various integrals, for example, Equation (6), a numerical 


integration was performed based on the Gauss-Legendre quadrature 
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method, whose operating equation is (Nieves 8 Domínguez, 1998; 
Campos-Aranda, 2003): 


A (19) 


2 


in which, w; are the method coefficients whose abscissas are h; and np the 
number of pairs in which the function f(x) is evaluated, with the argument 
indicated in f [-]. Davis and Polonsky (1972) obtained the 12 used pairs 
of w¡and h; with 15 digits, which are acceptable in the Basic language, as 
double precision variables. 


The operational equation of the bivariate integration is the following 
double integral (Nieves 8 Domínguez, 1998): 


fi 0 f(x, y) dxdy = A E A A E (20) 


Davis and Polonsky (1972) obtained the 20 pairs used in this double 
integral of w; and h; with 15 digits. 
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Association measures 
Concordance 


As the CF characterizes the dependence between the random variables u 
and v, it is necessary to study the measures of association, in order to 
have a method that allows estimating its parameter 6. In general terms, 
one random variable is concordant with another, when its large values are 
associated with the large values of the other and the small values of one 
with the small values of the other (Salvadori et al., 2007; Chowdhary € 
Singh, 2019). 


Some variables with direct linear correlation will be concordant, since 
when one increases the other also does. Variables with inverse linear 
correlation will be discordant, since large values of one will correspond to 
small values of the other. Then, (X1,y1) and (x2,y2) are said to be 
concordant if x¿<x2 and y1<y2, or x1>X2 and y1>y2. They are discordant 
ifxi<x2 and y1>y2, or x1>X2 and y1<y2. This implies that the pairs (x1- 
x2My1-y2)>0 are concordant (c) and discordant (d) when (x1-x2)M(y1- 
y2)<0 (Salvadori et al., 2007; Chowdhary € Singh, 2019). 


A numerical measure of association is a statistic that indicates the 
degree of dependence or association of the variables. For comparison 
purposes, such measures range from zero to +1 or -1, indicating the 
perfect positive association at +1 or negative at -1. Kendall's tau quotient 


and Spearman's rho coefficient are two nonparametric measures that 
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provide information on a special form of association or dependence, 


known as concordance (Salvadori et al., 2007; Chen € Guo, 2019). 


Kendall's tau quotient 


It measures the probability of having concordance pairs, so it is the 
quotient of c-d to c+d. The expression to estimate it with bivariate data 
is (Zhang € Singh, 2006; Zhang €: Singh, 2019): 


tn = E Ec signo — 3) 9) pena 


in the above equation, n is the number of observations and the sign[:] is 


+1 if such pairs are concordant and -1 if they are discordant. 


Genest and Favre (2007) present a test for the tau quotient, 
accepting that the null hypothesis Ho indicates that X and Y are 
independent and then the statistic has an approximately Normal 
distribution with zero mean and variance 2(2n+5)/[9n(n-1)]. Then, Ho 


will be rejected with a confidence level a = 5 % if: 


9n(n-1) 


Iznl > Za/2 = 1.96 (22) 


2(2n+5) 
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Spearman's rho coefficient 


It measures the correlation between pairs of ranks (R;,S;) of the random 
variables X; and Y;¡. Therefore, it is equivalent to the Pearson correlation 
coefficient. Its expression to estimate it in a bivariate record of size n is 
as follows (Chowdhary 8 Singh, 2019; Zhang € Singh, 2019): 


= Ri  B— (23) 


Pn = n(n+1D)(n-1) Hi=1 n-1 


Salvadori et al. (2007) present another formula for its calculation 
from a bivariate record with n data; which coincides with that used by the 
WMO (1971): 


n (p_o32 
_ 6 (Ri-Si) (24) 


n3—n 


Pn =1 


Genest and Favre (2007) also present a similar test for the rho 
coefficient, whose distribution is close to Normal with zero mean and 


variance 1/(n-1); therefore, Ho is rejected with a level a = 5 %, if: 


Vn —1- lp, 1 > 1.96 (25) 
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Estimation of the dependency parameter 


The simplest method to estimate the CF parameter 6 ¡is similar to the 
method of moments and is based on inverting the equation that relates O 
to Kendall's tau quotient or Spearman's rho coefficient. The other two 
available procedures are called: (1) maximum pseudo-likelinood method 
and exact maximum likelinood method (Meylan et a/., 2012; Chowdhary 
8: Singh, 2019; Zhang € Singh, 2019; Chen 8 Guo, 2019). To obtain 6, in 
equations (5), (10) and (13) we proceed by trial and error; on the other 


hand, in equations (3), (8) and (12) its value is cleared. 


Estimation of joint empirical probabilities 


The bivariate empirical probabilities were estimated based on the 
Gringorten formula, applied by Yue (2000b), Zhang and Singh (2019), 
and Chen and Guo (2019). The formula ¡s: 


E i-0.44 (26) 


Where ¡ is the number of each data, when they are ordered 
progressively and n is the total number of them, or amplitude of the 
processed register. The previous expression was applied in the two- 
dimensional plane, with the data ordered progressively; the maximum 


flows (Q) in the rows and the volumes (V) in the columns. The plane 
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formed is a square of n by n cells, with n cells on its main diagonal, when 
the order number of the row is equal to that of the column. Then each 
pair of annual data (Q and V) is located in the two-dimensional plane and 
the box defined by the intersection of the row and column is identified 


with the number / that corresponds to the historical year drawn. 


When the n pairs of data are drawn, year 1 is located and a 
rectangular or square area of minor Q and V values is defined, in which 
the count of numbered cells within, is NM or combinations of minor Q 
and V. Once the n values of NM; have been calculated, the Gringorten 
graphical position formula is applied to estimate the joint or bivariate 


empirical probability: 


F(xy)=P(Q<qV<wn===2 a 


n+0.12 


Selection of the Copula Function 


A simple Copula function selection approach is based on the fit error 
statistics, by comparing the observed empirical probabilities (w?) with the 
calculated theoretical ones (wf) with the CF being tested. The indicators 
applied are the standard mean error (EME, by its acronym in Spanish), 
the absolute mean error (EMA) and the maximum absolute error (EAM); 


their expressions are (Chowdhary € Singh, 2019): 
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EME = |=y7 ,(w? —w£)' (28) 


EMA = Yi, 1(w? — w£)| (29) 


EAM = max |(w? — wj)| (30) 


Dependency at the top end 


Generalities 


The most important criterion applied to select a CF is based on the 
dependency's magnitude at the upper end of the joint distribution, which 
has an impact on the extreme predictions” veracity. The upper right tail 
dependency (Au) is the conditional probability that Y is greater than a 
certain percentile (s) of FA(y), given that X is greater than that percentile 
in FxXQd); as s approaches the unit. The lower left tail dependence (A), 
compares that Y is less than X, when s approaches zero (Chowdhary 8 
Singh, 2019; Salvadori et al., 2007). 


In relation to the exposed CFs, those of Frank and Plackett, have 
insignificant dependencies in their extreme zones: therefore, A = O and 
Au = 0. In contrast, the Gumbel-Hougaard Copula has significant 


dependence on the upper tail, equal to: 
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Ay = 2 - 21/0 (31) 


Clayton and Raftery Copulas have it in their lower tail and equal to: 


dj 0 nta (32) 
.=2 (33) 


As already indicated, the Associated Clayton Copula has Ay = 271/80 
(Equation (17)). Families of unexposed copulas with Ay>0 are those of 
Galambos, Genest-Ghoudi and Húsler-Reiss (Salvadori et al., 2007; 
Chowdhary 8 Singh, 2019). 


Dupuis (2007) tested six Copula families and found that their ability 
to estimate extreme events ranges from poor to good, in the following 
order: Clayton, Frank, Normal, t-Student, Gumbel-Hougaard, and 
Associated Clayton (Survival Clayton). Similar conclusions are reached by 
Poulin et al. (2007), when comparing the same six families of Copulas and 


the one called A12, which has significant right tail dependence. 


Estimation of observed dependency 


To address the upper tail dependence (Ay) estimation shown by the 


available data, the Empirical Copula must be defined first. Since the CF 
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characterizes the dependence between the random variables X and Y, 
then the pair of ranges R; and S; from these variables are the statistic that 
retains the greatest amount of information and its scaling with the factor 
1/(n+1 ) generates a series of points in the unit square [0,1]?, forming 
the domain of the Empirical Copula (Chowdhary € Singh, 2019), defined 


as follows: 


C.(u,v) = 2311 (22 < 4, < v) (34) 


n+1 7 n+1 


In the above equation, 1(:*) denotes a function of the random 
variables U and V, which are a continuously increasing transformation of 
X and Y, relative to the empirical probability integrals F.(X) and FA(YN), 


whose equations are: 


__ Rango (X¡) 
TAR 


Rango (Y¡) 


Uj Ex) V === EY) (35) 


Poulin et al. (2007) use the estimator proposed by Frahm, Junker 
and Schmidt (2005), which is based on a random sample obtained from 


the Empirical Copula, its expression ¡s: 


1 1 1 1 
> =2- Zexp GXz in] inn An Gm (36) 
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This estimator accepts that the CF can be approximated by one of 
the extreme values class, it has the advantage of not requiring a threshold 
value for its estimation, like the four exposed by AghaKouchak, Sellars 
and Sorooshian (2013). 


Families of Copulas to be contrasted 
Justification 


Since the initial applications of Copula Functions (CF) in multivariate 
hydrological analyses, Favre et al. (2004) formulated in their Conclusions, 
that the crucial stage of such analysis lies in the adequate selection of the 
CFto the available data. Subsequently, Requena et al. (2013), highlighted 
that a CF that shows a good fit to the total available data, does not 
guarantee an adequate representation of the extreme values of such data 
and therefore, it is necessary to estimate the observed dependence A$“ 


(Equation (36)), for look for a CF that reproduces it. 


Therefore, when estimates of high joint return periods are required, 
greater than three times the amplitude of the available joint record of 
flows and maximum annual volumes, the CF that best reproduces the 
dependency on the right tail of the data (1%) must be sought and 
selected. In this study, four families of Copulas are contrasted, whose 
dependence on their right end varies from greater to lesser degree, these 


are: Joe, Associated Clayton, Gumbel-Hougaard and Student's t. The first 
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three belong to the class of Archimedean Copulas and the last one is of 
the elliptic class. 


Joe's Copula 


The basic equation of this CF family, which only accepts positive 
dependence, is the following (Joe, 1993; Salvadori et al., 2007; 
Chowdhary 8: Singh, 2019): 


C(u,v) =1-[U4-02+(41-w2-(1-w%4a-m.] 4 (37) 


where the dependency parameter is O > 1, with 0 = 1, for the case of 
independence between u and v. It does not have an expression that 
relates its parameter 0 to Kendall's tau quotient, therefore, it is estimated 
based on its generator equation g(s) and their derivative p'(s) of the CF, 


whose expression is (Salvadori et a/l., 2007; Chowdhary €: Singh, 2019): 


1 pts) 
Tn =1+4f4, rie (38) 


in which, s is the unitary random variable 0<s<1. The equations of its 


generator and its derivative are (Salvadori et a/., 2007; Zhang €: Singh, 
2019): 
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p(s) = —In[1 — (1 — sy?] (39) 
591 
¿=> (40) 


Taking the value of Kendall's tau quotient as data, Equation (38) 
was numerically integrated based on Equation (19), to obtain by trial and 
error the value of O that satisfies it. Obtained the value of 6, the significant 
dependence on the right end (Ay) of the Joe's Copula, is obtained by 


means of Equation (31) of the Gumbel-Hougaard Copula. 


Associated Clayton's Copula 


The Copula defined with Equation (18) will be applied and contrasted, 
since Dupuis (2007) exposed it and showed that it has a strong significant 
right end dependence (Au). 


Student't Copula 


Demarta and McNeil (2005) indicate that the Gumbel-Hougaard and 
Associated Clayton Copulas are simpler to apply than the extreme copulas 
derived from Student's t, but the latter defines a lower value of the 


dependency, as previously pointed out by Dupuis (2007) and confirmed 
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by Poulin et al. (2007). The Student's t-Copula equation is as follows 
(Salvadori et al., 2007; Chowdhary 8 Singh, 2019): 


dx - dy (41) 


pe 


_ ty (u), ty Mw) 1 x?-28-x-y+y? 
AC e 


In the above expression v > 2, O varies between -1 and 1, with 0 = 0 
for independence between u and v. t;*(:) is the inverse of the univariate 
Student's t-distribution with v degrees of freedom, for the probability 


indicated in parentheses. 


Then, with the idea of limiting the dependency in the right tail to a 
lower value, the Student t-Copula equations are exposed (Demarta 4 
McNeil, 2005; Poulin et a/., 2007; Salvadori et al., 2007; Genest, Favre, 
Beliveau, 8. Jacques, 2007): 


ay - arcseng (42) 
Ay = 2ty41 PEA (43) 


The minus sign inside the parentheses indicates that such value ¡s 
the exceedance probability (survival probability), which must define the 
random variable with Student's t distribution and v+1 degrees of 


freedom. To avoid the use of numerical algorithms that estimate the value 
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Of t,,,, based on the probability of exceedance, Demarta and McNeil (2005) 


present the following alternative formula: 


710/2 
2 = hn/2-arcseno)/a 05545 do 


1/2 
e cosVs:ds 


In the above equation, Á does not have a subscript, because the 
Student's t-Copula has the same significant dependence on its two tails 
or extremes. Equation (44) was numerically integrated based on Equation 
(19) and with these results, Table 1 was prepared, which is an extension 
of the one presented by Demarta and McNeil (2005). 


Table 1. Tail dependence values of the Student's t-Copula, for various 


values of v and 0. 


Association parameter value (0) 


-0.40 | -0.50 | -0.75 | 0.00 0.40 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 0.98 


0.0773 |0.0577| 0.0195 | 0.1817 | 0.3393 | 0.3910 |0.4195| 0.4502 | 0.4833 | 0.5195 | 0.5594 | 0.6042 | 0.6557 | 0.7177 | 0.7995 | 0.8729 


0.0378 | 0.0257 | 0.0061 | 0.1161 | 0.2606 | 0.3125 | 0.3419 | 0.3739 | 0.4091 | 0.4481 | 0.4918 | 0.5415 | 0.5995 | 0.6702 | 0.7648 | 0.8505 


0.0189 | 0.0117 | 0.0020 | 0.0756 | 0.2031 | 0.2532 | 0.2822 | 0.3144 | 0.3503 | 0.3907 | 0.4366 | 0.4896 | 0.5523 | 0.6298 | 0.7349 | 0.8311 


0.0049 | 0.0025 | 0.0002 | 0.0331 | 0.1269 | 0.1705 | 0.1970 | 0.2275 | 0.2625 | 0.3031 | 0.3506 | 0.4071 | 0.4758 | 0.5630 | 0.6845 | 0.7980 


0.0013 | 0.0006 | 0.0000 | 0.0150 | 0.0811 | 0.1173 | 0.1405 | 0.1679 | 0.2004 | 0.2393 | 0.2861 | 0.3434 | 0.4151 | 0.5086 | 0.6424 | 0.7699 


0.0004 | 0.0001 | 0.0000 | 0.0069 | 0.0527 | 0.0819|0.1015| 0.1255 | 0.1549 | 0.1911 | 0.2360 | 0.2925 | 0.3652 | 0.4627 | 0.6059 | 0.7452 


v = degrees of freedom. 
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Table 1 highlights that even for negative values and zero correlation, 
the Student's t-Copula has asymptotic dependence at ¡its extremes and 


that such dependence decreases with increasing degrees of freedom. 


The application or integration of the Student's t-Copula (Equation 
(41)) was carried out based on Equation (20), with 20 pairs of coefficients 
w; and abscissa hi; of 15 digits taken from Davis and Polonsky (1972) for 
the estimation of the integrals' upper limits, the algorithm exposed by 
Zelen and Severo (1972) was applied to estimate the random variable 
Student's t with v degrees of freedom, based on the Normal variable xp; 


in which, p is the probability of exceedance: 


£ xy + 2ebo) Solo) y galo) galo)... (45) 
being: 

gs) =0+x) (46) 
ga(x) = (5x5 + 1617 + 3x) (47) 
ga(x) = Gx? + 1915 + 17x? — 15x) (48) 
ga(x) = —(79x9 + 776x7 + 1482x5 — 1920x3 — 945x) (49) 


92160 
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When p < 0.50 make t = -t. The algorithm for estimating the 


Normal variable xp also comes from Zelen and Severo (1972) and is as 


follows: 


Co+c:8+c,:Ó? 
Xp = 9+C1 2 


in which: 


being: 

Co = 2.515517 
c1 = 0.802853 
cz =0.010328 
dí = 1.432788 
dz = 0.189269 
d3 =0.001308 


_ 1+d,:8+d2:82+d3:03 


(50) 


(51) 


If p > 0.50, make p =1 -— p and Xp = -Xp. The algorithms defined by 


equations (45) to (51) were tested with the numerical values exposed by 


Zhang and Singh (2019). 


To establish the 


integrals” lower limits of Equation (41), the 


tabulation of the Student's t distribution was used, as indicated by Ostle 
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and Mensing (1975) as Appendix 5. In said Table, the following values of 
t with v degrees of freedom from 2 to 10: 31.598, 12.941, 8.610, 6.859, 
5.959, 5.405, 5.041, 4.781 and 4.587. Such values with a minus sign 
were taken as lower limits according to the degree of freedom assigned 
in the application of Equation (41); which, logically, varied from 2 to 10; 
selecting the Student's t-Copula of best fit (equations (28) to (30)). 


Ratification of the selected Copula function 


This is the most important stage in the process of the CF practical 
application, since it is verified that the model faithfully reproduces the 
observed joint probabilities (Equation (27)). Yue (2000a) indicates that 
the simple and practical way of representing the empirical and theoretical 
joint probabilities consists of taking the first one to the abscissa axis and 
the second to the ordinate axis; logically, in such a graph, each pair of 
data defines a point that coincides with or moves away from the line at 
450, The inspection of the graph described and the value of the correlation 
coefficient, in these cases, greater than 0.98, ratify the validity of the 


joint probabilistic model. 


Yue (2000b), and Yue and Rasmussen (2002) apply the Kolmogorov- 
Smirnov test with a significance level (a) of 5 %, to accept or reject the 
maximum absolute difference (dma) between the joint probabilities. To 
evaluate the statistic (D,) of the test, the expression proposed by Meylan 
et al. (2012), for a = 5 % this is: 
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D, = A (52) 


n is the number of data. If the dma is less than D,, the adopted CF is 
ratified. 


Bivariate frequency analysis concepts 


Bivariate return period 


The first bivariate return period of the event (X, Y) is defined under the 
OR condition, which indicates that the limits x ó y, or both can be 
exceeded and then, the classical equation of the return period or inverse 


probability of exceedance will be (Shiau et a/., 2006; Genest 8: Chebana, 
2017): 


1 1 1 


E ara) TN - EARIRON E 


where C[Fx(0), Fy(y)] is the selected CF. 


The second bivariate return period of the event (X, Y) is associated 
with the case in which both limits are exceeded (X > x, Y > y) or AND 


condition, its equation is (Shiau et al., 2006; Genest €: Chebana, 2017): 
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PASAS) Fx y Goy) 1-FxG0)-FyO)+C[FxGO,Fy()] 


Aldama (2000) obtains the expression F;xy(x,y) of the bivariate 
probability of exceedance through a simple and logical reasoning of 
probabilities applied in the Cartesian plane. Instead, Yue and Rasmussen 
(2002) use the Cartesian plane to define the bivariate event (X, Y) 


conceptually, which can occur in any of the four quadrants. 


The relationship between bivariate and univariate return periods is 
as follows (Yue 8 Rasmussen, 2002; Shiau et al., 2006; Vogel 8 
Castellarin, 2017): 


Txy < min[7y, Ty] < max[Tx, Ty] < T'xy (55) 
being: 
1 1 
Xx“ FG) 1-Fx() (38) 


_ 1 _— 1 
Y EO)  1-Fy0) 


(57) 


Salvadori and De Michele (2004) introduce in detail the concept of 


the Secondary Return Period (7), designated in this way to emphasize that 
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the joint return period Txy is the primary one, from where it comes when 


using the ¡solines defined by the applied CF, whose expression is: 
L.=|u,w er ce) =8 (58) 


where s is the unitary random variable 0<s<1 and C is the CF tested. 
Then, a region Bc(s) is defined in the unit space (1?) above the i¡soline, 


below it and to the left, which will be: 
Bels) = [(u, v) € 1?:C(u,v) < s) (59) 


In a CF of the Archimedean class, the Kendall univariate 
distribution, designated Kc(s), provides a measure of the events within 
the Bc(s); ¡its equation is (Salvadori € De Michele, 2004; Salvadori 8. De 
Michele, 2007; Salvadori et al., 2007; Gráler et al., 2013): 


Ke(s) =s-— + (60) 


in which, g(s) is the generator of the CF and g'(s) its derivative. Finally, 


the secondary return period (7) of the events outside of Bc(s), is: 


¿=== (61) 
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where the denominator is the probability of exceedance (survival 
function), which corresponds to possibly of destructive or dangerous 
events. The relationship between (€ and the joint return periods is as 
follows (Gráler et a/., 2013): 


a (62) 


The application of Equation (61) requires the prior estimation of 
Equation (60), for the following six non-exceedance probabilities (s): 
0.98, 0.99, 0.998, 0.999, 0.9998 and 0.9999, which correspond to the 
joint return periods of 50, 100, 500, 1 000, 5 000 and 10 000 years. In 
the Gumbel-Hougaard and Joe Copulas, the equations of its generator 


and its derivative are the following (Zhang and Singh, 2019): 


p(s) = Cinsy? (63) 
9) =ÍEmya (64) 
p(s) = —In[1 — (1 — sy?] (65) 
p'(5) = 2 (66) 
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Critical or design events 


Volpi and Fiori (2012) highlight that the plot of the AND-—type joint return 
period, shown later as Figure 2, presents a severe inconsistency as it 
contains, in a bivariate context, univariate critical thresholds. Due to the 
above, such a graph is considered to be made up of two portions, the two 
designated simple (naive part) and the correct one (proper part). The 
straight parts are the tails or straight lines asymptotes to the curved part. 
The probability of an event occurrence or pair of Q and V, is variable in 
the curved part and decreases along the straight part, although all the 
values define the same joint return period. In summary, the pairs of 
values of the asymptote lines have low probabilities of occurrence and 
therefore should not be included in the analysis of the search for critical 
or severe floods (Q and V). For practical purposes, the extreme points of 
the curved part can be defined according to their empirical distribution or 


close to the beginning of the asymptotic lines. 


Selection of marginal distributions 


The approach for the selection of the marginal distributions was very 
simple and consisted of applying the three PDFs that have been 
established as reference or applicable under precept, which are the Log- 
Pearson type Ill (LP3), the General Extreme Values (GVE) and 


Generalized Logistics (LOG). In addition, three widely used distributions 


Tecnología y ciencias del agua, ISSN 2007-2422, 


Open Access bajo la licencia ee BY-NC-SA 4.0 14(5), 120-188. DOI: 10.24850/j-tyca-14-05-03 
(https: //creativecommons.org/licenses/by-nc-sa/4.0/) 


a 1) Check for updates 
OPEN ACCESS 


Tecnología y 


CienciaszAgua 


were applied: the Generalized Pareto (PAG), the Kappa and the Wakeby. 
The first four mentioned PDFs have three fit parameters and the last two, 
four and five. With the exception of LP3, which was applied with the 
moments method, in the logarithmic (WRC, 1977) and real (Bobée, 1975) 
domains, the rest were fitted with the L-moments method (Hosking 8 
Wallis, 1997). 


The selection of the most convenient FDP was based on the value 
that each one generates with a non-exceedance probability of 1 %; that 
is, a very low value that does not exceed the minimum values observed, 
to avoid negative marginal probabilities (u,v). This selection criterion is 
the one required in records of Q and V that present very low minimum 


values, compared to their maximum extremes. 


In addition, the standard errors of fit (EEA) and mean absolute errors 
(EAM) were taken into account, as well as the predictions” magnitude in 


return periods greater than 500 years. 


Fitting errors 


The first criteria applied for the selection of the best PDF to some available 
data or series were the fitting errors (Kite, 1977; Willmott € Matsuura, 
2005; Chai 8, Draxler, 2014). This criterion and the one described to avoid 
negative probabilities allow an adequate distribution to be adopted 
between the models: LP3, GVE, LOG, PAG, Kappa and Wakeby. 


Changing the probabilities observed by the ordered data of the 


analyzed series (x; or y;) in equations (28) and (29), and the probabilities 
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calculated by the estimated values with the PDF that is tested or 
contrasted, the standard error of fit (EEA) and the mean absolute error 
(EAM) are obtained. The estimated values (2, or y,) are searched for the 
same probability of non-exceedance, assigned to the data by Gringorten's 


empirical formula (Equation (26)). 


Data to be processed 


In order to compare the four CfFs that have a significant right tail 
dependence (Ay), the 16 joint records of maximum flow (Q) and volume 
(V) of the annual floods reported by Aldama et al. (2006), which add up 
to 620 pairs of data, in 15 important reservoirs in Mexico and one in 


project (La Parota). 


In addition, the joint record of annual Q and V of the floods recorded 
at the La Cuña hydrometric station, exposed by Gomez et al. (2010) and 
Campos-Aranda (2022). The record of Q and V of the floods at the 
entrance to the Malpaso Dam, shown by Domínguez and Arganis (2012), 


was also processed. 


Table 2 shows the record of the 26 joint data of Q and V of the floods 
at the entrance to the Adolfo Lopez Mateos Dam on the Humaya River in 
the state of Sinaloa, Mexico. Such record comes from Aldama et al. (2006) 


and processed as a numerical application. 
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Table 2. Maximum flows, volumes and their bivariate order numbers of 
the 26 annual floods at the entrance to the Adolfo López Mateos Dam 
(Humaya), Mexico (Aldama et al., 2006). 


Año E A NM; | Año e nl NM; 
(m3/s) | (Mm?) (m3/s) | (Mm) 
1974 600 140 4 | 1987 401 160 6 
1975 290 200 4 | 1988 305 150 3 
1976 316 70 3 | 1989 581 170 8 
(1977 | 242 | so | 1 | 1990 | 2035 | 1150 | 22 | 
(1978 248 | 420 | 4 | 1991| 488 | 190 | 9 | 
1979 | 891 | 400 | 15 | 1992 | 242 | 310 | 3 | 
(1980 461 | 190 | 7 | 1993| 487 | 160 | 7 | 
(1981 | 9245 | 4580 | 26 | 1994 | 933 | 150 | 6 | 
¡1982 | 2388 | 1410 | 23 | 1995 | 127 | 70 | 1 | 
(1983 | 277 | 160 | 3 | 1996 | 4490 | 1240 | 24 | 
1984 481 290 9 | 1997 904 160 9 
1985 614 450 16 | 1998 | 2529 1020 | 22 
1986 1064 550 20 | 1999 | 1193 380 17 
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Wald-Wolfowitz Test 


This nonparametric test has been used by Bobée and Ashkar (1991), Rao 
and Hamed (2000), and Meylan et al. (2012) to verify independence and 
stationarity in records of maximum annual flows (X;). Therefore, it was 


proposed to apply this test to the joint records of Q and V, which must be 
random samples. 


Results and their discussion 
Contrast of Copula functions with 2,>0 


Randomness verification 


Table 3 shows, in progressive order of basin size (A in km), the 16 joint 
records of Q and V processed, taken from Aldama et al. (2006) and at the 
end, the two joint records also processed. Column 4 shows the number of 
data (n) of each record and columns 5 and 6, the value of the U statistic 
of the Wald-Wolfowitz Test; whose absolute value must be less than 1.96 
for such registration to be random. 
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Table 3. General data and results of the Wald-Wolfowitz Test in the 18 


joint records of Q and V processed. 


A WW Test 
No. Name of the dam or record n 
(km?) Q V 
1 Madin 171 30 | -0.757 | -0.335 
2 Eustaquio Buelna (Guamúchil) 1630 37 | (4.156) | (2.597) 


3 Josefa Ortiz de Domínguez (El Sabino) 2 268 34 | (1.986) | 1.529 


4 Abelardo L. Rodríguez 2 430 44 | (2.220) | (2.346) 
5 Sanalona 3 250 38 | 0.337 | 1,555 
6 Ignacio Allende (La Begoña) 4 984 27 | 0.287 1.241 
7 La Parota 7 067 31 | -0.760 | -0.532 
8 Benito Juárez (El Marqués) 9 697 26 | -0.601 0.427 
9 Adolfo Ruiz Cortines (Mocúzari) 10 719 61 0.147 1.189 
10 Adolfo López Mateos (Humaya) 10 972 26 0.677 1.252 


11 Belisario Domínguez (La Angostura) 18 537 33 0.955 0.247 


12 La Boquilla 21 003 64 | -0,.084 | -0.251 
13 Luis Donaldo Colosio (Huites) 26 020 52 | -0.091 | -0.602 
14 Venustiano Carranza (Don Martín) 31 034 52 0.726 1.054 
15 Plutarco Elías Calles (El Novillo) 58 280 22 | 0.637 0.496 
16 Adolfo López Mateos (Infiernillo) 108000 | 43 | -0.707 | -0.022 
17 Estación hidrométrica La Cuña 19097 | 55 | 0.284 0.213 
18 Presa Malpaso (Netzahualcóyotl) 34 800 47 | -0.666 | 0.558 
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It can be deduced from the cited values that two joint records are 
not random, those of the Guamúchil and Abelardo L. Rodríguez dams. 
Specific statistical tests were applied to such records and it was found in 
the Guamúchil dam that the Q and V series show persistence and excess 
variability due to their extreme values. In addition, the Q series shows an 
upward trend. In the Abelardo L. Rodriguez dam, both series of Q and V 
show persistence, excess variability and a negative trend. The Q record 
of the El Sabino dam can be accepted randomly, since it only shows slight 
persistence. Due to the above, the joint records of the Guamúchil and 
Abelardo L. Rodríguez dams were eliminated from the subsequent 
analyses. 


Values of association measures 


In column 2 of Table 4, the Kendall's tau quotient values obtained with 
Equation (21) are displayed. Regarding Spearman's rho coefficient, shown 
in column 2 of Table 5, the values that do not have parentheses were 
equal to the calculated with equations (23) and (24). Instead, the 
quantities in parentheses come from Equation (23), which leads to slightly 
larger values. 
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Table 4. Kendall's tau quotient values of the observed dependency 
(A5F“) and of the applied Joe and Clayton Copulas Associated with the 


indicated joint record. 


Dam or Joe Copula Clayton Asoc. 
Ta A A A 
record 0 dy 0 y 
| Madin  [|0.6598| 0.7356 | 4.687 | 0.8406 | 3.8784 | 0.8363 | 


El Sabino 0.5936 | 0.6776 | 3.748 | 0.7969 | 2.9211 0.7888 


Sanalona 0.5733 | 0.6053 | 3.525 | 0.7827 | 2.6867 0.7726 


La Begoña 0.6866 | 0.7831 | 5.177 | 0.8567 | 4.3818 | (0.8537) 


La Parota 0.5226 | 0.6724 | 3.039 | 0.7438 | 2.1892 | (0.7286) 


El Marqués [0.5385 | 0.6462 | 3.179 | 0.7564 | 2.3333 0.7430 


Mocúzari 0.7093 | 0.8024 | 5.630 | 0.8690 | 4.8797 | (0.8676) 


Humaya 0.4523 | 0.6834 | 2.522 | 0.6837 | 1.6517 0.6573 


La Angostura | 0.4318 | 0.4730 | 2.397 | 0.6647 | 1.5200 0.6338 


La Boquilla 0.6091 | 0.6440 | 3.938 | 0.8075 | 3.1168 | 0.8006 


Huites 0.4495 | 0.5655 | 2.504 | 0.6811 | 1.6329 0.6541 


Don Martín 0.6456 | 0.7446 | 4.455 | 0.8317 | 3.6426 0.8267 


El Novillo 0.4286 | 0.6166 | 2.378 | 0.6616 | 1.5000 0.6300 


Infiernillo 0.4064 | 0.4714 | 2.253 | 0.6398 | 1.3694 0.6028 


La Cuña 0.7199 | 0.7819 | 5.887 | 0.8750 | 5.1394 | 0.8738 


Presa Malpaso | 0.3747 | 0.4159 | 2.092 | 0.6072 | 1.1982 0.5608 
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Table 5. Spearman's rho coefficient values of the observed dependency 
(A5F“) and of the Gumbel-Hougaard Copulas and Student's t in the 


indicated set of records. 


Dam or Gumbel-Hougaard t de Student 
Pn Ago 
record 0 w w(v=2) | w(v= 10) 
Madin (0.8603) | 0.7356 | 2.9392 0.7340 0.6678 0.3834 
El Sabino 0.6883 0.6776 | 2.4605 0.6746 0.6072 0.2965 
Sanalona 0.7595 0.6053 | 2.3433 0.6558 0.5889 0.2725 


La Begoña 0.8107 0.7831 | 3.1909 | (0.7574) | 0.6929 0.4225 


La Parota 0.6968 0.6724 | 2.0946 | (0.6077) | 0.5443 0.2184 


El Marqués 0.7648 0.6462 | 2.1667 0.6230 0.5581 0.2345 


Mocúzari (0.8827) | 0.8024 | 3.4399 | (0.7768) | 0.7142 0.4570 


Humaya (0.7258) | 0.6834 | 1.8258 | 0.5383 0.4849 0.1563 


La Angostura 0.5836 0.4730 | 1.7600 | 0.5173 0.4681 0.1409 


La Boquilla 0.7922 0.6440 | 2.5584 | 0.6888 0.6212 0.3157 


Huites 0.5980 0.5655 | 1.8164 | 0.5354 0.4825 0.1541 


Don Martín (0.8303) | 0.7446 | 2.8213 DE72I15 0.6546 0.3636 


El Novillo 0.5788 0.6166 | 1.7500 0.5140 0.4654 0.1385 


Infiernillo (0.5707) | 0.4714 | 1.6847 0.4910 0.4476 0.1233 


La Cuña (0.9088) | 0.7819 | 3.5697 0.7857 0.7243 0.4737 


Presa Malpaso | (0.5510) | 0.4159 | 1.5991 0.4574 0.4226 0.1038 
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Dependency contrast (A, ) 


Columns 3 of Table 4 and Table 5 show the values of the observed 
dependence (AF“), estimated with Equation (36). The CF selection 
criterion was that its significant dependence (Ay), should be somewhat 
larger or slightly smaller than the observed dependency (25%). With this 
criterion in Table 4, the Joe CFs were defined for the record of the Adolfo 
Lopez Mateos Dam (Humaya) and the Associated Clayton CF for the 
Plutarco Elias Calles Dam (El Novillo) record. 


On the other hand, in Table 5, with the selection criteria described, 
10 records are defined whose dependence on their right tail can be 
presented by the Gumbel-Hougaard CF (GH) and three that are included 
between the Associated Clayton CF and GH, whose dependency values 
(Au) are shown in parentheses. The CF of the Student's t reproduces the 


dependence of records from the Malpaso Dam, with v= 2. 


Search of the marginal distributions 
Distribution of maximum annual flow 


Table 6 shows the fit errors and predictions (m3/s) obtained with the six 
distributions applied to the record of maximum flows in Table 2, whose 
minimum value recorded is 127 m/s. Of these six PDFs, the only ones 


that define values with a probability of non—-exceedance of 1 % that are 
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lower than the aforementioned flow are the GVE and the LOG. Taking into 


account the fit errors, the GVE is adopted. 


Table 6. Fit errors and predictions (m3/s) of the six distributions applied 
in the record of maximum annual flows of the floods entering the Adolfo 


Lopez Mateos Dam (Humaya), Mexico. 


Return periods in years 
EEA EAM 


PDF 
(m3*/s) | (m?/s) 50 100 500 1 000 5 000 10 000 


LP3 230.8 132.6 9 062 15 192 48 746 79778 | 247 228 | 400 939 


GVE 675.5 279.8 6 663 10 381 28 653 44 244 | 121000 | 186 452 


LOG 745.8 301.9 6 364 9 961 28 041 43 766 123060 | 192 052 
PAG 631.9 255.6 6 824 10 412 26 828 40 002 100 278 | 148 626 
KAP 585.1 243.0 7072 10 697 26 739 39 253 94 582 137 656 


WAK 676.7 286.5 6 770 10 379 27 111 40 694 103672 | 154 793 


The location (us), scale (ai) and shape (k1) parameters of the 
adopted GVE distribution are: 447.7865, 372.8246, -0.6236228, where 


the equation ¡s: 


F(x) = exp (- [1 — ki cp) (67) 


a1 
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Distribution of annual volumes 


In Table 7, similar to the previous one for the annual volumes of Table 2, 
whose minimum value is 50 Mm, Again, the only PDFs that do not 
generate negative marginal probabilities are the GVE and LOG. The GVE 


is adopted because it shows less fit errors. 


Table 7. Fit errors and predictions (m3/s) of the six distributions applied 
to the record of annual volumes of the floods at the entrance to the 


Adolfo López Dam Mateos (Humaya), Mexico. 


EEA EAM Return periods in years 
ds (m?/s) | (m*/s) 50 100 500 1 000 5 000 10 000 
LP3 | 1898 | 96.6 | 4061 | 6726 | 20745 | 33252 97 487 153 958 
GVE | 370.6 | 138.8 | 3021 | 4726 | 13163 | 20401 56 269 87 019 
LOG | 400.4 | 144.6 | 2881 | 4528 | 12856 | 20134 57 053 89 331 
PAG | 355.1 | 133.1 | 3091 | 4739 | 12345 | 18 494 46 869 69 796 
KAP | 338.6 | 135.1 | 3207 | 4875 | 12327 | 18186 44 336 64 877 
WAK | 375.3 | 146.4 | 3075 | 4729 | 12428 | 18697 47 873 71 630 


The location (uz), scale (a2) and shape (k2) parameters of the 
adopted GVE distribution are: 193.4795, 167.2692 and -0.6288097, 


where the expression is: 
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m0) =exp(-[1-k,222]%) (68) 


Selection and ratification of the CF 


The bivariate data processing in Table 2 led to the following three 
indicators of association: rxy = 0.9670, 1, = 0.4523 and Pn = 0.6738. 
Equation (22), with n = 26 and the cited tau, gives a value of 3.24; 
therefore, the Kendall quotient is significant. Equation (25) generates a 


value of 2.26, so Spearman's coefficient is also significant. 


On the other hand, Table 8 shows the statistical fit indicators that 
were obtained by applying the Clayton, Frank, Plackett, t-Student, 
Gumbel-Hougaard (GH), associated Clayton and Joe CFs. In equations 
(28) to (30), the empirical bivariate probabilities were estimated with 
Equation (27) and the theoretical ones with equations (2), (4), (9), (41), 
(7), (18) and (37). The values reported for the Copula t-Student 


correspond to a v = 6, where the fit errors have already stabilized with 


0 = sen (==) = 0.6522, according to Equation (42). 
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Table 8. Statistical indicators of the fit of the Copula functions indicated 


in the annual input floods to the Adolfo López Mateos Dam (Humaya), 


Mexico. 

Copula 0 EME EAM | No. DP | No. DN | MDP MDN Ay 

Clayton(C.) | 1.6517 | 0.0529 | 0.0365 14 12 0.1505 | -0.0779 | 0.0000 
Frank 4.9310 0.0461 | 0.0321 15 dd 0.1334 | -0.0701 | 0.0000 
Plackett 11.079 | 0.0441 | 0.0315 11 15 0.1291 | -0.0727 | 0.0000 
t-Student 0.6522 | 0.0482 | 0.0347 13 3 0.1415 | -0.0643 | 0.2641 
GH 1.8258 | 0.0473 | 0.0354 14 1,2 0.1352 | -0.0601 | 0.5383 
C. Asociada | 1.6517 | 0.0475 | 0.0362 14 42 0.1277 | -0.0549 | 0.6573 
Joe 2.5220 | 0.0477 | 0.0365 AS 11 0.1267 | -0.0533 | 0.6837 


New acronyms: 


DP, DN = positive and negative differences. 


MDP, MDN = maximum positive and negative difference. 


Based on the results in Table 8, there is no difficulty in selecting the 


Plackett CF as the best CF for the data in Table 2, since it reports the 


smallest statistical indicators (EME and EAM) (shown in italics). 


However, when taking into account that the observed dependency 


(A5F9) reached a value of 0.6834, the CF that should be selected is Joe, as 
already stated in Table 4. Table 8 shows that the fit of Joe's CF to the data 
is better than the poorest of all, which was defined by Clayton's CF. It is 


noteworthy that Joe's CF gives rise to the lowest values of the maximum 
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differences between empirical and theoretical probabilities (shown in 


italics and shaded). 


Table 9 shows the bivariate probabilities of non-exceedance, 


observed empirical (w?) and calculated theoretical (w?) with Joe's CF. The 


maximum positive and negative differences are also shaded. 


Table 9. Joint non—-exceedance probabilities and their differences, 


calculated with Joe's CF, for the floods annual entrance to the Dam 


Adolfo Lopez Mateos (Humaya), Mexico. 


No. wi wi Differences | No. wi wj Differences 
1 | 0.1363 | 0.1896 0.0533 14 0.2129 0.1619 0.0509 
2 |0.1363 | 0.1300 0.0063 15 | 0.0980 | 0.1051 -0.0071 
3 |0.0980 | 0.0311 0.0669 16 | 0.2894 | 0.2392 0.0502 
4 |0.0214 | 0.0101 0.0114 17 | 0.8254 | 0.8643 0.0389 
5 |0.1363 | 0.1378 -0.0015 18 | 0.3277 | 0.2369 0.0909 
6 | 0.5574 | 0.5663 0.0089 19 | 0.0980 | 0.1215 0.0234 
7 (0.2511 | 0.2249 0.0263 20 | 0.2511 | 0.1955 0.0556 
8 |0.9786 | 0.9851 0.0066 21 | 0.2129 | 0.2460 -0.0331 
9 |0.8637 | 0.8939 0.0302 22 | 0.0214 | 0.0051 0.0163 
10 | 0.0980 | 0.0973 0.0007 23 | 0.9020 | 0.9195 -0.0175 
11 | 0.3277 | 0.3210 0.0067 24 | 0.3277 | 0.2666 0.0612 
12 | 0.5957 | 0.4690 0.1267 25 | 0.8254 | 0.8767 0.0513 
13 | 0.7489 | 0.6690 0.0798 26 | 0.6340 | 0.6058 0.0282 
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On the other hand, Equation (52) defines D, = 0.2663 and since the 
maximum absolute difference in Table 9 is 0.1267, the Kolmogorov- 
Smirnov test ratifies the adopted Joe's CF. The correlation coefficient (rxy) 
between the empirical and theoretical probabilities, estimated with Joe's 
CF, was 0.9897; therefore acceptable. The graphic contrast between both 
probabilities, to ratify its adoption, is shown in Figure 1 for the data in 
Table 9. 


órica 


r 
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0.0 0.2 0.4 0.6 0.8 1.0 
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Figure 1. Graphical contrast of joint probabilities (Q, V) calculated with 
Joe FC de Joe, for annual floods at the entrance of the Adolfo Lopez 


Mateos (Humaya) dam, Mexico. 
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Bivariate return period graphs Ty, 


The bivariate return periods of the AND type are estimated with Equation 
(54). For their estimations with values of Tyy 0f 500, 1 000, 5 000 and 10 
000 years, Joe's CF was applied. Volumes and peak flows are arbitrarily 
selected to obtain their marginal (equations (67) and (68)) and joint 
(Equation (54)) probabilities of non-exceedance. Table 10 shows results 


to define the four graphs of Figure 2. 


Table 10. Couples of peak flow and annual volume used to define the 
graphs of the joint return period type AND with Joe's CF, in the floods 
input to the Adolfo López Mateos Dam (Humaya), Mexico. 


Tr Ti Tr Tr 

500 years 1 000 years 5 000 years 10 000 years 

Vv Q v Q v Q v Q 
Mm? m3/s Mm? m3/s Mm? m3/s Mm? m3/s 

0 28 653 0 44 244 0 121 000 O | 186452 | 
1000 | 28636 | 2000 | 44200 | 4000 | 120940 | 5000 | 186420 | 
2000 | 28571 | 4000 |44020 | 8000 | 120735 | 10000 | 186 300 | 
3000 | 28439 | 6000 | 43642 | 12000 | 120290 | 15000 | 185 840 | 
5000 | 27903 | 8000 | 43002 | 16000 | 119500 | 20000 | 185 250 
6000 | 27453 | 10000 42014 | 20000 | 118390 | 25000 | 184 150 
7000 | 26836 | 12000 | 40542 | 24000 | 116830 | 30000 | 182 800 


Tecnología y ciencias del agua, ISSN 2007-2422, 
14(5), 120-188. DOI: 10.24850/j-tyca-14-05-03 


4.0 


Open Access bajo la licencia cc BY-NC-SA 
(https: //creativecommons.org/licenses/by-nc-sa/4.0/) 


OPEN Que: 


0) Check for updates 


Tecnología y 
CienciaszAgua 
Tor Tr Tor Tor 
500 years 1 000 years 5 000 years 10 000 years 
v Q v Q v Q v Q 
Mm? m3/s Mm? m3/s Mm? m3/s Mm? m3/s 

8000 | 26000 | 13000 | 39560 | 28000 | 114700 | 35000 | 181 120 | 
9000 | 24858 | 14000 | 38 348 | 32000 | 111900 | 40000 | 178750 | 
10000 | 23266 | 15000 | 36850 | 36000 | 108 140 | 45000 | 175800 | 

10500 | 22222 | 16000 | 34960 | 40000 | 103060 | 50000 | 172 050 

11000 | 20940 | 17000 | 32520 | 42000 | 99850 | 52500 | 169 960 

11500 | 19315 | 18000 | 29212 | 44000 | 96060 | 55000 | 167 520 

12000 | 17165 18500 | 27050 | 46000 | 91430 |57500 | 164 740 

12 500 | 14046 | 19000 | 24335 | 47000 | 88770 |60000 | 161 570 

12750 | 11745 19500 | 20720 | 48000 | 85800 |62500 | 157 970 

13000 | 8125 | 19750 | 18330 | 49000 | 82460 | 65000 | 153 850 

13050 | 6990 | 20000 | 15170 | 50000 | 78670 |67500 | 149 310 

13100 | 5490 | 20100 | 13530 | 51000 | 74250 |70000 | 143 810 

13 163 0 20200 | 11520 | 52000 | 69140 | 72000 | 138 820 

20300 | 8730 | 53000 | 62840 |74000 | 133 200 

20350 | 6450 | 53500 | 59250 |76000 | 126 530 

20 401 0 54 000 | 54985 | 78000 | 118 720 

54 500 | 50160 |80000 | 109 000 

55 000 | 44100 | 81000 | 103 580 

55 500 | 35840 |82000 | 96950 

56000 | 23600 |83000 | 89650 
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Tx Txy Txy Txy 
500 years 1 000 years 5 000 years 10 000 years 
v Q v Q v Q v Q 
Mm3 m3/s Mm m3/s Mm m3/s Mm m3/s 
56 269 0 84 000 80 650 
85 000 68 990 
86 000 53 500 
87 019 0 
100 
90-7'(Q ,V)=10,000 años 
An” 80 
€ 
= 70 
m 
O  60-T'(Q ,V)=5,000 años 
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7 50 
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ce 40 
1) 
S 30 
T'(Q ,V )=1,000 años 
E 20 (Q V) 
o 
> 10 
0 
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Figure 2. Graphs of the four design joint return periods Ty, obtained 
with Joe's CF, in the annual floods at the entrance to the Adolfo Lopez 


Dam Mateos (Humaya), Mexico. 
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Selection of design events 


In Figure 2 or in Table 10, infinite pairs of Q and V can be selected, which 
satisfy the joint design return period and which are defined as a subgroup 
of critical pairs, as they are within the curved portion of each 7;, graph, 
outside the asymptote lines (Volpi € Fiori, 2012). 


The combinations of peak flow and volume that have the same design 
bivariate return period, establish floods or hydrographs that will produce 
different effects in the reservoir that is designed or revised; adopting for 
security, the one that generates the most critical conditions. To form each 
design hydrograph, there are theoretical and empirical methods (Aldama, 
2000; Aldama et al., 2006; Campos-Aranda, 2008; Gómez et al., 2010; 
Gráler et al., 2013). 


Contrast of the Joint Return Periods 


Table 11 and Table 12 show the calculations made to verify equations 
(55) and (62), based on the Gumbel-Hougaard and Joe CFs. It is observed 
and verified in all cases that both CFs comply with Equation (55). For the 
secondary return period, Equation (62) is exclusively fulfilled with Joe's 


CF in the last three return periods analyzed. 
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Table 11. Estimated univariate and bivariate return periods with the 
Gumbel-Hougaard CF in the floods at the entrance to the Adolfo Lopez 


Mateos Dam (Humaya), Mexico. 


Gumbel-Hougaard CF Txy secondary 
Tx y Ty 
CIFx(),Fr(y)] Txy Txy Kc(s) y 

50 0.9709003 34.4 91.7 0.990844 109.2 

100 0.9854162 68.6 184.6 0.995450 219.8 

500 0.9970778 342.2 927.9 0.999094 | 1 104.2 
1 000 0.9985386 684.3 1856.7 | 0.999547 | 2 209.6 
5 000 0.9997076 3 420.4| 9294,.9 | 0.999910 | 11 059.5 
10 000 0.9998538 6 836.7 |18 600.0 | 0.999955 | 22 104.4 


Table 12. Estimated univariate and bivariate return periods with the Joe 


CF in the floods at the entrance to the Adolfo Lopez Mateos Dam 


(Humaya), Mexico. 


Tx y Tv Joe CF Txy secondary 
CIFX(O,Fr(y) 1 Txy Txy Kc(s) E 
50 0.9736739 38.0 23:L 0.987932 82.9 
100 0.9868368 76.0 146.3 0.993994 166.5 
500 0.9973674 379.9 731.4 0.999022 | 1 022.7 
1 000 0.9986837 759.7 1 462.7 | 0.999000 | 1 000.0 
5 000 0.9997368 3799.2| 7313.5 — 0.999800 | 5 000.7 
10 000 0.9998683 7 594.9 | 14 627.0 | 0.999900 | 9 998.3 
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Conclusions 


In Mexico, the six works cited: Aldama (2000), Ramírez-Orozco and 
Aldama (2000), Escalante-Sandoval and Reyes-Chávez (2002), Aldama 
et al. (2006), Gómez et al. (2010), and Domínguez and Arganis (2012), 
highlight the efforts made to develop reliable flood estimation methods 
that provide security to rural and urban hydraulic works that are 
dimensioned and/or revised with them. When seeking to estimate 
bivariate floods, as approximate as possible; this study adds to the 


aforementioned efforts. 


The fundamental advantage of using the Copula Functions (CF) in 
bivariate floods frequency analysis (AFCb, by its acronym in Spanish), lies 
in easily constructing the joint probability distribution, based on the same 
or different marginal univariate distributions adopted; prior estimation of 
the dependence between the random variables: maximum flow (Q) and 


runoff volume (V) from the available annual floods. 


The application of the CF in the AFCb has shown that the selection 
of the ideal Copula is based on reproducing the dependence of the 
extreme right observed; which guarantees a more approximate estimate 
of the predictions associated with high joint return periods. Therefore, CFs 
that exhibit, from plus to minus, various degrees of right tail dependency, 
are now tested by the families of Joe's Copulas, Associated Clayton, 


Gumbel-Hougaard, and Student's t. 


The contrast made with 18 annual records of Q and V collected, 


detected two non-random ones, which were eliminated. For the rest, Joe's 
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CF is suggested in a record, processed as a numerical application. The 
Associated Clayton CF is adopted in another registry and the Gumbel- 
Hougaard CF is proposed in 10 registries. In three registers, its observed 
dependency is between the one exhibited by these two mentioned CfFs. 
Finally, it was found that the CF of the Student's t reproduces the 


dependence of the Malpaso Dam record, with v=2 (see Table 5). 


The CFs that do not show such dependence, such as the Clayton, 
Frank and Plackett, are also applied to the available Q and V data and 
allow validating the quality of the fit achieved by the tested CF. The 
foregoing is supported by the fact that all the CFs cited and exposed show 


a similar fit to the data (see Table 8). 


The processing of the 26 annual data of Q and V of the floods at the 
entrance to the Adolfo Lopez Mateos Dam (Humaya), in the state of 
Sinaloa, Mexico; It allowed to show the numerical development of the 
exposed theory and to describe in detail Figure 1 and Figure 2, which 
include the acceptance of the CF adopted and the basic graph of the 


results or predictions of Q and V of the joint return period type AND. 
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